function fun = define_function_fun(par)

%%%% We define functions here with parameters given by par
%%%% We assume log-normal distribution for productivity as in the data

    %%%% Log-normal PDF and CDF; notice that the lower bound is close to zero (instead of 0) 
    %%%% and the upper bound is very large (instead of infinity) were used to speed up computation.
    fun.f  = @(x, a) 1 ./ (x .* par.log_eps_sig * sqrt(2 * 3.14159265359))...
        .* exp(- (log(x) - par.log_eps_mu - a).^2 / 2 / par.log_eps_sig.^2)...
        .* (x >= par.EPS_L_BOUND) .* (x <= par.EPS_H_BOUND);
    fun.F  = @(x, a) logncdf(x, par.log_eps_mu + a, par.log_eps_sig) .* (x >= par.EPS_L_BOUND) .* (x <= par.EPS_H_BOUND);
    fun.l_bound     = @(a)  par.EPS_L_BOUND ;
    fun.h_bound     = @(a)  par.EPS_H_BOUND ;

    %%%% Buyer threshold productivity given seller productivity x; 
    %%%% check eq. (14) in the paper
    %%%% K_K_tilde is defined as K / K_tilde; L is a small part of Psi_0
    %%%% and see the result of IID case in eq. (23).
    fun.threshold    = @(x, L, CHI_PI, CHI_PI_k, K_K_tilde) L ./ (1 - par.THETA -  CHI_PI - CHI_PI_k * K_K_tilde)...
        - par.THETA ./ (1 - par.THETA -  CHI_PI - CHI_PI_k * K_K_tilde) .* x; 

    %%%% The following two functions assume x as buyer productivity, y as seller productivity
    %%%% Note: the order of seller and buyer is VERY important!
    
    %%%% Unit surplus for buyer (after considering credti). The denominator of eq. (15)
    fun.Delta        = @(y, x, B, CHI_q, CHI_PI) max(0.000001, (1 - par.tau_k) .* B .* ((1 - par.THETA - CHI_PI) .* x + par.THETA .* y) + (1 - par.DELTA) .* (1 - CHI_q));
    
    %%%% Marginal gain of one unit of capital (see eq (18))
    fun.capital_gain = @(y, x, a_sell, a_buy) (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy) ; 

    %%%% Marginal gain of one unit of money (see eq (18))
    fun.money_gain   = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI) (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        ./ fun.Delta(y, x, B, CHI_q, CHI_PI); % marginal gain of one unit of money

    %%%% Marginal profit gain of one unit of money (see eq (18))
    fun.profit_gain   = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI) (1 - par.tau_k) .* B .* x .* (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        ./ fun.Delta(y, x, B, CHI_q, CHI_PI); 

    %%%% Liquidity gain (in reallocating capital in partial sale transactions)
    fun.liq_gain   = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K) (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        .* (Z_K + CHI_k .* (1 - par.DELTA) + CHI_PI_k .* (1 - par.tau_k) .* B .* x) ./ fun.Delta(y, x, B, CHI_q, CHI_PI);

    fun.reall_L      = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k)...
        (fun.Delta(y, x, B, CHI_q, CHI_PI) + CHI_q .* (1 - par.DELTA) + CHI_PI .* (1 - par.tau_k) .* B .* x) .* fun.f(y, a_sell) .* fun.f(x, a_buy);

    fun.reall_P      = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K)...
        (fun.Delta(y, x, B, CHI_q, CHI_PI) + CHI_q .* (1 - par.DELTA) + CHI_PI .* (1 - par.tau_k) .* B .* x) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        .* (Z_K + CHI_k .* (1 - par.DELTA)  + CHI_PI_k .* (1 - par.tau_k) .* B .* x) ./ fun.Delta(y, x, B, CHI_q, CHI_PI);

    fun.reall        = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, K_tilde_K)...
        (fun.Delta(y, x, B, CHI_q, CHI_PI) + CHI_q .* (1 - par.DELTA) + CHI_PI .* (1 - par.tau_k) .* B .* x) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        .* min((Z_K + CHI_k .* (1 - par.DELTA) + CHI_PI_k .* (1 - par.tau_k) .* B .* x) ./ fun.Delta(y, x, B, CHI_q, CHI_PI), K_tilde_K);

    % fun.surplus_buy =  @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, L, K_tilde_K) (1 - par.tau_k) .* B .* (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
    %     .* min((L .* (1 - par.tau_k) .* B + (1 - par.tau_k) .* CHI_PI .* B .* x + (1 - par.DELTA) .* (1 - CHI_q)) ./ fun.Delta(y, x, B, CHI_q, CHI_PI), K_tilde_K);

    fun.surplus_buy =  @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, K_tilde_K) (1 - par.tau_k) .* B .* (x - y) .* fun.f(y, a_sell) .* fun.f(x, a_buy)...
        .* min((Z_K + CHI_k .* (1 - par.DELTA) + CHI_PI_k .* (1 - par.tau_k) .* B .* x) ./ fun.Delta(y, x, B, CHI_q, CHI_PI), K_tilde_K);


    % Quantity in transaction
    fun.q_fun    = @(y, x, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde)...
        min(1, (Z_K_tilde + (1 - par.DELTA) .* CHI_k + CHI_PI .* (1 - par.tau_k) .* B .* x) ./ max( (1 - par.tau_k) .* B .* ((1 - par.THETA - CHI_PI).* x + par.THETA .* y)...
       + (1 - par.DELTA) .* (1 - CHI_q), 1e-6));

    % using cash first in meetings; does not change allocation
    % fun.debt_fun = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_k, Z_K) fun.f(y, a_sell) .* fun.f(x, a_buy)...
    %     .* max(0,  fun.q_fun(y, x, B, CHI_q, CHI_PI, CHI_k, Z_K) .* ( (1 - par.tau_k) .* B .* ((1 - par.THETA - CHI_PI) .* x + par.THETA .* y) + (1 - par.DELTA) ) - Z_K);

    % using debt first in meetings
    fun.debt_fun = @(y, x, a_sell, a_buy, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde) fun.f(y, a_sell) .* fun.f(x, a_buy)...
        .* min(...
        fun.q_fun(y, x, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde) .* ( (1 - par.tau_k) .* B .* ((1 - par.THETA - CHI_PI) .* x + par.THETA .* y) + (1 - par.DELTA) ),...
        fun.q_fun(y, x, B, CHI_q, CHI_PI, CHI_k, Z_K_tilde) .* ( (1 - par.DELTA) * CHI_q + (1 - par.tau_k) .* B .* CHI_PI) + (1 - par.DELTA) * CHI_k + (1 - par.tau_k) .* B .* x .* CHI_PI...
        );

    %%%% For computing productivity dispersion
    fun.M = @(x, L, a)fun.F(x, a) - integral(@(eps_t)...
        (fun.F(max(L / (1 - par.THETA) -  par.THETA / (1 - par.THETA) * eps_t, eps_t), a) - fun.F(eps_t, a)) .* fun.f(eps_t, a),...
        fun.l_bound(a), x);
    fun.m = @(x, L, a) fun.f(x, a) - fun.f(x, a) .* (fun.F(max(L / (1 - par.THETA) -  par.THETA / (1 - par.THETA) * x, x), a) - fun.F(x, a));
    fun.G = @(x, L, a) fun.M(x, L, a) / fun.M(fun.h_bound(a), L, a);
    fun.g = @(x, L, a) fun.m(x, L, a) / fun.M(fun.h_bound(a), L, a);

end